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The  numerical  tools  typically  used  to  model  the  evolution  of  fluid  instabilities  in  inertial  confinement  fusion 
(ICF)  hydrodynamics  codes  are  examined,  and  some  are  found  to  have  properties  which  would  seem  to  be 
incompatible  with  the  accurate  modeling  of  small-amplitude  perturbations,  i.e.,  perturbations  in  the  linear  stage 
of  evolution.  In  particular  a  “differentiability  condition”  which  is  satisfied  by  the  physics  in  such  situations 
is  not  necessarily  satisfied  by  the  numerical  algorithms  in  typical  use.  It  is  demonstrated  that  it  is  possible  to 
remove  much  of  the  non-differentiability  in  many  cases,  and  that  substantial  improvement  in  one’s  ability  to 
accurately  model  the  evolution  of  small  amplitude  perturbations  can  result.  First  a  simple  example  involving  a 
non-differentiable  radiation  transport  algorithm  is  shown,  and  then  the  non-differentiabilities  introduced  by  the 
use  of  upwind  and  “high  resolution”  hydrodynamics  algorithms  are  analyzed. 


I.  INTRODUCTION:  NUMERICALLY  MODELING  THE 
EVOLUTION  OF  SMALL  AMPLITUDE  PERTURBATIONS 

In  Fig.  1  we  show  a  schematic  diagram  of  the  problem 
which  motivates  this  paper.  We  seek  to  numerically  model 
the  laser-driven  implosion  of  a  nearly  spherically  symmetric 
layered  target  in  such  a  way  as  to  effect  the  thermonuclear 
burn  of  of  its  fusible  components.  In  particular,  we  wish  to  be 
able  to  model  the  evolution  of  perturbations  away  from  spher¬ 
ical  symmetry  with  sufficient  accuracy  to  be  able  to  reliably 
compute  the  expected  energy  gain  of  the  pellet.  These  per¬ 
turbations,  which  can  grow  via  Rayleigh-Taylor,  Richtmyer- 
Meshkov,  and  other  fluid  instabilities,  can  originate  in  the  ini¬ 
tial  conditions  of  the  pellet,  or  be  impressed  upon  the  pellet  by 
inhomogeneities  in  the  laser  radiation  field.  If  they  are  small 
enough,  they  will  initially  undergo  a  period  of  linear  evolu¬ 
tion.  If  they  are  large  enough,  they  may  eventually  transition 
to  a  more  and  more  nonlinear  regime,  perhaps  culminating  in 
driven,  fully  turbulent  flow.  We  are  interested  here  in  the  accu¬ 
rate  modeling  of  the  early,  linear,  stages  of  perturbation  evo¬ 
lution,  for  this  evolution  forms  the  critical  initial  conditions 
inherited  by  the  nonlinear  regime. 

Although  the  inertial  confinement  fusion  (ICF)  pellet  im¬ 
plosion  motivates  our  discussion,  we  are  actually  interested 
in  a  much  more  general  question:  What  are  the  appropriate 
numerical  tools  to  accurately  model  the  evolution  of  small- 
amplitude  perturbations,  when  the  unperturbed  state  itself 
may  be  time-dependent,  highly  nonlinear,  and  may  contain 
spatial  and  temporal  scales  that  are  smaller  than  we  can  actu¬ 
ally  resolve,  e.g.,  shock  waves? 

It  is  important  to  note  that  in  the  discussion  throughout 
this  paper  we  are  assuming  that  all  of  the  physical  features 
of  interest  are  smooth  on  some  spatial  and  temporal  scales, 
albeit  scales  that  may  be  smaller  than  we  can  afford  to  actu¬ 
ally  resolve  numerically.  Thus  when  we  pose  questions  of  the 
physics,  we  assume  differentiability  at  some  finite  scale  ev- 
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FIG.  1:  Geometry  of  a  laser-driven  inertial  confinement  fusion  (ICF) 
pellet.  We  are  interested  in  accurately  modeling  the  evolution  of  per¬ 
turbations  to  this  nominally  spherically-symmetric  configuration. 


erywhere,  but  when  we  pose  numerical  questions,  we  do  not 
necessarily  assume  that  all  such  scales  can  be  resolved  on  any 
mesh  that  we  can  afford  to  use.  Even  if  we  are  treating  inter¬ 
faces  between  materials,  we  are  assuming  either  that  there  is 
smoothness  at  some  very  small  scale,  or,  equivalently,  that  the 
imposition  of  smoothness  on  some  small  scale  will  not  affect 
the  answers  to  any  questions  in  which  we  have  interest. 

Since  the  above  assumption  of  smoothness  at  some  finite 
scale  may  cause  some  readers  to  pause,  lets  us  examine  the 
reasons  that  we  make  such  an  assumption  here,  as  well  as 
the  plausibility  of  the  assumption.  The  reader  will  see  be¬ 
low  that  we  predicate  our  analysis  on  the  linear  behavior  of 
additive  infinitesimal  perturbations,  both  in  in  the  initial  con¬ 
ditions  and  in  the  subsequent  evolution  of  the  system.  But  a 
system  that  contains  interfaces,  i.e.,  true  discontinuities  in  the 
state  variables,  either  in  the  initial  conditions  or  in  its  subse¬ 
quent  evolution,  simply  cannot  be  described  in  terms  of  ad- 
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ditive  perturbations  unless  one  is  using  a  coordinate  system 
that  is  moving  with  the  interfaces,  e.g.,  a  Lagrangian  coordi¬ 
nate  system  in  the  case  of  material  interfaces.  Thus  it  could 
be  argued  that  the  analysis  we  present  below  is  not  applica¬ 
ble  to  such  situations  unless  one  is  using  such  a  coordinate 
system.  Nonetheless  we  believe  strongly  that  the  analysis  we 
present  below  is  applicable.  Our  reasons  are  two.  First,  there 
are  physical  dissipative  processes  active  in  all  plasmas  which 
prevent  the  maintenance  of  discontinuities  in  state  variables, 
e.g.,  mass  diffusion,  viscosity,  and  conductivity.  The  scale 
sizes  produced  by  these  mechanisms  may  be  smaller  than  can 
be  resolved  numerically,  in  which  case  a  numerical  simula¬ 
tion  that  included  these  mechanisms  would  be  virtually  iden¬ 
tical  to  one  which  did  not,  but  the  physical  scale  sizes  are 
finite  nonetheless,  thus  making  our  assumptions  valid.  Sec¬ 
ond,  the  interface  instabilities  in  which  we  are  interested  in 
general  do  not  have  their  behavior  significantly  changed  if  we 
approximate  them  as  narrow  but  smooth  transitions  between 
two  states.  This  is  explicitly  acknowledged  in  the  formulae  we 
use  to  approximate  the  growth  rate  of  the  ablative  Rayleigh- 
Taylor  instability,  for  example,  wherein  a  factor  (1  -|-  kL)  1 
pre -multiplies  the  Rayleigh-Taylor  interface  growth  rate  term 
to  allow  for  finite  plasma  gradient  scale  length  effects,  where 
k  is  the  transverse  Fourier  wave  number  of  the  perturbation 
and  L  is  the  gradient  scale  length  ([1],  pg.  62).  Clearly  a  suf¬ 
ficiently  small  but  finite  L  has  a  negligible  effect  on  the  inter¬ 
face  growth  rate  in  this  instance.  Another  example  is  the  fact 
that  accounting  for  physical  viscosity  and  heat  conduction  in 
the  compressible  fluid  equations  near  a  shock  wave  broadens 
the  shock  front,  but  does  not  change  the  shock’s  propagation 
speed  or  the  jumps  across  the  shock  front.  (This  is  even  true 
numerically  as  well  as  physically  if  we  use  numerical  meth¬ 
ods  which  conserve  mass,  momentum  and  total  energy  at  the 
discrete  level,  as  we  do  here.)  Nonetheless,  the  reader  should 
be  cognizant  of  the  fact  the  the  arguments  we  present  below 
cannot  be  rigorously  shown  to  be  applicable  to  situations  in 
which  true  discontinuities  exist. 


II.  THE  LINEAR  EVOLUTION  OF  INFINITESIMAL 
PERTURBATIONS:  THE  DIFFERENTIABILITY 
CONDITION 

We  will  frame  our  discussion  within  that  of  systems  of  con¬ 
servation  laws,  which  take  the  form 

+y.f,-(ffi ,q2,...,qm,x,t)  =  0;  i=l,m  (1) 

where  x  and  t  are  space  and  time  respectively,  and  m  is  the 
number  of  conservation  laws  in  the  system.  Examples  of 
such  equations  include  the  Navier-Stokes  equations,  the  equa¬ 
tions  of  magnetohydrodynamics  (MHD),  the  Vlasov  equation, 
passively-driven  convection,  and,  of  course,  our  ICF  pellet  im¬ 
plosion  scenario. 

As  we  have  stated  before,  we  shall  assume  that  the  initial 
conditions  q(x.O)  are  smooth  at  some  sufficiently  small  spa¬ 
tial  scale,  and  that  there  is  present  in  the  physics  some  dissipa¬ 
tive  mechanism  active  at  some  finite  but  perhaps  small  spatial 


and  temporal  scales  to  ensure  that  the  solution  q(x,t)  is  also 
smooth  at  those  scales.  We  do  not,  however,  assume  that  we 
can  actually  resolve  all  such  scales  numerically. 

The  above  notwithstanding,  we  will  try  to  impose  on  our 
numerical  algorithms  the  following  property  of  the  physics: 
If  we  confine  our  perturbation  to  that  of  the  initial  condition 
t  =  0,  then  in  the  limit  of  vanishing  perturbation  amplitudes 
about  virtually  any  unperturbed  time-evolving  system,  there 
exists  a  constant  of  proportionality  between  any  component 
of  the  initial  perturbation  and  any  other  component  of  the 
evolved  perturbation  at  some  fixed  later  time.  This  constant 
of  proportionality  is  precisely  the  partial  derivative  of  the  lat¬ 
ter  with  respect  to  the  former. 

In  particular  the  following  “linearity”  property  holds:  We 
define  a  physical  time  evolution  operator  T  which  advances 
the  physical  solution  in  time  from  t  =  0  to  t  =  tf. 

q{tf)  =  T(q(  0))  (2) 

where  we  have  dropped  the  notational  dependence  of  q  on  x 
for  clarity.  We  do  the  same  for  the  perturbation  quantities  p  \ 
and  p2  below. 

For  the  laser-driven  ICF  problem  of  interest  here,  T  is  a 
nonlinear  differentiable  operator  on  q ,  i.e.,  q(t f)  is  a  differen¬ 
tiable  function  of  q(0),  unless  we  find  ourselves  at  a  bifurca¬ 
tion  point  in  chaotic  flow.  For  the  purposes  of  this  paper  we 
assume  that  we  are  not  at  such  a  bifurcation  point. 

We  now  introduce  infinitesimal  perturbations  p  \  (x)  and 
P2(x)  to  the  initial  conditions.  From  the  differentiability  of 
T  and  the  assumption  that  p  \  (x)  and  p2  (x)  are  infinitesimal  it 
follows  that  in  the  limit  of  vanishing  perturbation  amplitudes: 

T(q+  api)  -T(q)  =  a(T(q  +  p{)  -T(q))  (3) 

where  a  is  a  scalar,  and 

T(q  +  Pl+P2)~T(q)  =  T(q  +  pl)  +  T(q  +  p2)-2T(q)  (4) 

where  we  have  dropped  the  notational  dependence  of  p  \  and 
P2  on  x  for  clarity. 

Equations  (3)  and  (4)  above  define  a  perturbation  evolution 
that  is  normally  described  as  “linear.”  They  apply  whether  the 
growth  or  decay  of  the  perturbation  is  exponential  in  time,  as 
it  is  for  the  Rayleigh-Taylor  instability,  or  has  some  other  time 
dependence,  as  is  the  case  for  the  Richtmyer-Meshkov  insta¬ 
bility.  The  terminology  aside,  they  describe  the  physical  be¬ 
havior  of  infinitesimal  perturbations  to  the  basic  flow  defined 
by  T(q).  It  is  the  primary  thesis  of  this  paper  that  a  neces¬ 
sary,  or  at  least  highly  desirable,  condition  for  the  success¬ 
ful  numerical  modeling  of  the  evolution  of  extremely  small 
perturbations  to  a  specified  physics  problem,  like  that  of  the 
imploding  ICF  pellet,  is  that  the  numerical  time  evolution  op¬ 
erator  T„  satisfy  the  very  same  equations.  Henceforth,  when  a 
given  numerical  time  evolution  operator  Tn  satisfies  equations 
(3)  and  (4),  we  shall  say  that  that  Tn  satisfies  the  “differentia¬ 
bility  condition.”  Surprisingly  perhaps,  many  of  the  numerical 
methods  in  use  in  modern  ICF  codes  fail  this  test. 


3 


III.  NON-DIFFERENTIABLE  ALGORITHMS  IN  MODERN 
ICF  CODES 

In  the  last  section  we  put  forth  the  thesis  that  a  numer¬ 
ical  time  evolution  operator  Tn  should  satisfy  the  differen¬ 
tiability  condition,  i.e.,  satisfy  equations  (3)  and  (4),  if  we 
wish  to  accurately  model  the  evolution  of  small  perturba¬ 
tions  to  a  physical  system.  A  legitimate  question  for  any¬ 
one  not  fully  familiar  with  ICF  codes  might  be  the  follow¬ 
ing:  How  could  such  non-differentiability  possibly  find  its 
way  into  such  codes?  Looking  at  Eq.  (1),  we  see  that  our 
numerical  task  is  simply  to  perform  discretizations  of  spatial 
and  temporal  derivatives.  Certainly  the  obvious  finite  differ¬ 
ence  and  finite  element  formulations  of  such  derivatives  do 
not  introduce  any  non-differentiability,  so  barring  any  inher¬ 
ent  non-differentiability  in  the  equations  themselves,  where 
could  such  non-differentiability  possibly  originate?  Rather 
than  attempting  to  answer  this  question  comprehensively,  let 
us  simply  list  some  of  the  places  where  one  might  find  a  lack 
of  differentiability  in  such  codes: 

1.  “High  resolution”  numerical  algorithms  for  fluid  dy¬ 
namics,  sometimes  known  as  “modern  front-capturing 
methods”  or  ‘modern  shock-capturing  methods.” 

2.  Upwind  methods  for  fluid  dynamics  (at  sonic  points) 

3.  “Sharp  cutoff"  flux  limiters  for  thermal  conductivity, 
and  for  radiation  diffusion 

4.  Linear  interpolation  in  table  lookups  for  equation  of 
state  data  and  opacities 

The  non-differentiability  associated  with  the  last  item  above 
will  hopefully  be  obvious  to  the  reader,  as  will  the  fact  that 
creating  a  table  lookup  interpolation  procedure  with  contin¬ 
uous  derivatives,  that  also  avoids  the  possibility  of  spurious 
interpolated  values,  may  be  a  challenging  problem  for  some 
data  sets.  The  remaining  three  items  will  be  explored  in  the 
following  sections. 

A.  Symptoms  of  Numerical  Non-differentiability 

An  infinitesimally  perturbed  solution  that  encounters  a 
point  of  non-differentiability  will  exhibit  a  false  nonlinear  re¬ 
sponse  of  some  kind.  A  pure  Fourier  mode  in  a  periodic 
problem  may  not  remain  pure,  for  example.  This  can  only 
happen  if  the  numerical  algorithm  used  has  points  of  non¬ 
differentiability  inside  an  “envelope”  of  the  solution  space 
spanned  by  the  unperturbed  solution  on  the  one  hand  and 
by  its  perturbed  counterpart  on  the  other.  In  that  case,  the 
perturbed  solution  must  at  some  point  in  its  history  have 
“crossed”  a  point  of  non-differentiability,  perhaps  many  times 
if  the  perturbation  is  oscillatory  in  space  or  time,  seeing  one 
derivative  on  one  side  of  the  crossing  point,  and  another  on  the 
other.  Clearly  as  the  amplitude  of  the  perturbation  becomes 
smaller  and  smaller,  the  envelope  in  solution  space  in  which 
the  point  of  non-differentiability  must  reside  becomes  smaller 


and  smaller.  Thus,  unless  a  point  of  non-differentiability  re¬ 
sides  directly  on  the  unperturbed  solution,  one  can  easily  be¬ 
come  just  plain  “lucky”  in  any  given  calculation,  with  no 
symptoms  at  all,  and  yet  encounter  a  major  “glitch”  with  the 
same  numerical  software  on  a  modestly  different  problem.  We 
will  see  an  example  of  this  in  a  later  section.  On  the  other 
hand,  even  if  a  point  of  non-differentiability  is  encountered, 
the  result  may  or  may  not  be  harmful,  or  even  measurable. 
This  will  clearly  depend  on  many  factors,  including  the  size 
and  location  of  the  point  of  non-differentiability. 

Thus,  by  simply  running  computational  examples  with  var¬ 
ious  algorithms,  as  we  do  here,  it  will  be  impossible  to 
prove  that  the  reason  an  algorithm  failed  is  that  it  was  non- 
differentiable,  or  that  the  reason  another  succeeded  was  that  it 
was  differentiable.  All  we  can  do  here  is  try  to  lend  some  em¬ 
pirical  plausibility  to  our  hopefully  plausible  hypothesis.  Our 
task  is  made  all  the  more  difficult  by  the  fact  that  it  is  prob¬ 
ably  impractical  to  construct  real  ICF  codes  with  algorithms 
that  are  completely  differentiable  over  all  of  solution  space. 

B.  An  Example  of  the  Effects  of  Numerical 
Non-differentiability 

As  we  have  stated  previously,  one  of  the  possible  sources 
of  non-differentiability  in  ICF  codes  is  the  use  of  sharp  cut¬ 
off  radiation  diffusion  and  thermal  conductivity  flux  limiters. 
We  refer  the  reader  to  ([2],  pp  478-481)  for  the  motivation  and 
form  of  such  flux  limiters  in  the  radiation  diffusion  case,  the 
case  we  shall  examine  in  this  section.  Briefly,  the  diffusion 
approximation  used  in  many  ICF  codes  can  produce  fluxes  of 
radiant  energy  F(s)  that  exceed  what  is  physically  possible 
Fmax(s).  Here  s  represents  all  of  the  parameters  on  which  F 
may  depend,  e.g.,  the  local  radiation  energy  density  gradient 
and  Rosseland  mean  opacity.  In  that  case,  one  of  the  remedies 
is  to  limit  the  flux  F(s)  so  that  it  does  not  exceed  Fmax(s).  The 
obvious  “sharp  cut-off"  solution  is  simply 

E(s)  =  rnin(E(s),Fm“'(s))  (5) 

Clearly  Eq.(5)  defines  a  flux  which  is  a  non-differentiable 
function  of  its  arguments  at  all  values  of  s  for  which  F(s)  = 
Fmax(s).  An  alternative  flux  limiter  is  the  so-called  “harmonic 
mean”  limiter: 


F(s)  =  F(s)/(l+F(s)/Fmax(s))  (6) 

which  is  a  differentiable  function  of  its  arguments.  These  are 
the  two  flux  limiters  that  we  examine  below. 

In  Figures  2  and  3,  we  show  the  results  of  a  pair  of  2D 
r  —  9  laser  pellet  implosion  simulations  using  the  Naval  Re¬ 
search  Laboratory’s  (NRL’s)  FAST  code  [3].  The  pellet  and 
imposed  laser  intensity  profile  are  described  in  [4].  Briefly, 
the  pellet  is  a  pure  deuterium-tritium  (DT)  shell  of  outer  ra¬ 
dius  0.169  cm  and  thickness  0.034  cm,  driven  by  0.351  jlm 
wavelength  laser  radiation  at  an  intensity  of  10  TW  for  the 
first  4  nsec,  and  ramping  monotonically  to  450  TW  at  6.05 
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FIG.  2:  Legendre  mode  amplitude  of  areal  mass  density  vs.  time  for 
a  FAST  run  of  a  NIF  pellet  using  an  extremely  small  Legendre  1  =  2 
perturbation,  and  using  a  harmonic  mean  (differentiable)  radiation 
flux  limiter.  The  Legendre  spectrum  should,  and  is,  dominated  by 
the  nearly  static  seed  perturbation  during  this  early-time  evolution. 
The  seed  mode  is  plotted  as  a  thick  line  with  embedded  “+”  symbols, 
while  the  other  Legendre  modes  are  plotted  as  simply  thin  lines 


nsec,  where  it  stays  until  laser  turn-off  at  8.6  nsec.  A  single 
1  =  2  Legendre  mode  perturbation  of  extremely  small  ampli¬ 
tude  is  imposed  on  the  outer  surface.  The  quantity  of  interest 
is  the  integrated  (in  r)  areal  mass  density  as  a  function  of  0, 
here  analyzed  in  Legendre  space.  We  calculate  only  the  very 
early  time  behavior  here,  terminating  the  simulation  shortly 
after  the  first  main  shock  has  broken  out  through  the  interior 
surface  of  the  pellet.  The  physics  is  such  that  one  should  see 
a  nearly  static  evolution  of  the  1  =  2  mode,  with  little  or  no 
generation  of  other  Legendre  modes  during  this  early  phase 
of  the  pellet  implosion.  In  Fig.  2  we  show  Legendre  mode 
amplitude  vs.  time  for  a  calculation  which  used  the  harmonic 
mean  radiation  diffusion  flux  limiter  Eq.  (6)  ,  which  is  a  dif¬ 
ferentiable  function  of  its  arguments.  Note  that  the  behavior  of 
the  Legendre  modes  is  qualitatively  what  we  expect:  a  nearly 
static  seed  mode,  with  small  amplitudes  maintained  for  all  the 
other  modes.  In  Fig.  3  we  show  the  results  for  a  simulation 
which  was  identical  to  that  in  Fig.  2,  except  that  we  are  now 
using  the  sharp  cut-off  version  of  the  radiation  diffusion  flux 
limiter  Eq.  (5),  which  is  a  not  a  differentiable  function  of  its 
arguments.  Note  the  dramatic  difference  in  the  results,  with 
the  unphysical  generation  of  higher  order  harmonics  which 
eventually  drive  even  the  primary  l  =  2  mode  away  from  its 
physically  correct  amplitude. 

This  one  simple  example  makes  it  clear  that  non- 
differentiable  numerical  algorithms  are  a  potential  serious 
threat  to  the  accuracy  of  simulations  of  small-amplitude  per¬ 
turbations  on  ICF  laser  pellets. 


IV.  UPWIND  METHODS  AND  “HIGH  RESOLUTION” 
METHODS  AS  SOURCES  OF  NON-DIFFERENTIABILITY 

While  we  again  wish  to  emphasize  that  any  source  of  nu¬ 
merical  non-differentiability  should  be  viewed  as  a  threat  to 


FIG.  3:  Same  as  Fig.  2,  but  using  a  sharp  cut-off  (non-differentiable) 
radiation  flux  limiter.  Note  the  quick  non-physical  rise  of  other 
Legendre  modes,  a  numerical  artifact  caused  by  the  use  a  non- 
differentiable  numerical  algorithm. 


accuracy,  we  will  for  the  remainder  of  this  paper  focus  on 
hydrodynamics  algorithms.  Specifically,  we  will  investigate 
the  first  two  sources  of  non-differentiability  from  our  list  in 
the  previous  section:  “high  resolution”  methods  and  upwind 
methods.  The  use  of  these  methods  as  opposed  to  standard 
(differentiable)  finite  difference  and  finite  element  methods 
derives  from  the  same  need:  the  need  to  robustly  deal  with 
structures  too  small  to  be  resolved  on  any  reasonable  spatial 
or  temporal  mesh,  as  we  explore  in  the  paragraphs  that  follow. 

As  we  stated  in  our  assumptions,  as  long  as  there  is  some 
dissipative  mechanism  active  at  some  finite  but  perhaps  small 
scales,  solutions  to  (1)  with  smooth  initial  conditions  will  re¬ 
main  smooth  for  all  time.  However,  in  many  such  systems, 
including  compressible  fluid  dynamics,  the  system  can  and 
often  does  evolve  in  such  a  way  to  produce  small  regions  over 
which  the  solution  or  its  derivatives  can  change  dramatically. 
If  these  regions  are  smaller  than  a  mesh  size,  from  a  numer¬ 
ical  viewpoint  we  are  treating  discontinuities.  We  will  term 
such  unresolved  regions  “fronts”  for  the  purpose  of  this  paper. 
Given  that  we  do  not  wish  to  actually  resolve  these  fronts,  we 
demand  instead  that  they  be  represented  numerically  as  nar¬ 
row  regions  on  our  grid  without  numerical  artifacts  that  con¬ 
taminate  adjacent  regions,  that  they  propagate  with  the  correct 
speed,  and  that  the  jumps  across  them  are  physically  correct. 
In  most  cases  this  situation  is  addressed  by  the  Lax-Wendroff 
Theorem  ([5]),  which  states  that  if  one’s  numerical  approxi¬ 
mation  to  Eq.  (1)  is  in  “flux”  or  “conservation”  form,  and  the 
solution  converges  everywhere  but  on  a  set  of  measure  zero  to 
some  solution,  and  in  addition  satisfies  an  entropy  inequality, 
then  the  above  demands  will  be  met.  Thus  the  great  majority 
of  methods  designed  to  treat  fronts  in  the  context  of  Eq.  (1) 
are  in  conservation  form,  i.e.,  a  form  consisting  of  numerical 
fluxes  connecting  adjacent  grid  points,  these  fluxes  being  used 
to  advance  the  numerical  solution  in  time. 

Using  conservation  form  does  not  by  itself  give  the  desired 
result  however,  since  one  still  needs  to  compute  a  convergent 
solution.  In  general,  numerical  methods  not  designed  to  deal 
with  fronts  will  not  produce  the  desired  convergence  in  their 
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presence,  often  producing  numerical  oscillations  that  degrade 
the  solution  severely.  Most  often,  successful  methods  use 
some  sort  of  explicit  or  implicit  artificial  numerical  dissipa¬ 
tion  to  achieve  such  convergence. 

Over  the  past  30  years  a  host  of  algorithms,  known  vari¬ 
ously  as  “modern  front-capturing  methods,”  “modern  shock¬ 
capturing  methods,”  or  “high  resolution  methods,”  have  been 
developed  in  an  attempt  to  perform  calculations  containing 
fronts  more  accurately  and  more  efficiently  than  with  the  more 
traditional  artificial  dissipation  approaches.  The  first  of  these 
methods  was  flux-corrected  transport  (FCT)  [6,  7],  but  there 
are  now  a  large  number  of  others  (e.g,  MUSCL  [8],  PLM  [9], 
PPM  [10],  ENO  [11],  WENO  [12,  13],  TVD  [14]  methods, 
and  Discontinuous  Galerkin  (DG)  methods  [15]  ).  What  dis¬ 
tinguishes  the  “high  resolution”  methods  from  their  predeces¬ 
sors  is  their  attempt  to  constrain  the  numerical  fluxes,  grid 
point  by  grid  point  and  timestep  by  timestep,  in  such  a  way 
as  to  avoid  the  production  of  unphysical  values  in  the  solution 
vector  q. 

A  well-known  theorem  by  Godunov  [16]  states  that  any 
linear  algorithm  that  guarantees  the  avoidance  of  unphysical 
values  renders  the  algorithm  at  most  first  order  accurate  in 
time  and  space.  Thus  all  modern  front-capturing  algorithms 
of  order  greater  than  one  must  be  inherently  nonlinear  opera¬ 
tors.  However,  almost  without  exception,  these  modern  front¬ 
capturing  methods  involve  the  use  functions  that  are  not  just 
nonlinear,  but  non-differentiable  as  well,  using  min,  max,  abs, 
and  sign  operators,  as  well  as  similar  functions  involving 
if-then-else  statements.  Thus  at  a  fundamental  level,  time 
evolution  operators  incorporating  such  algorithms  are  non- 
differentiable  functions  of  their  arguments.  Compounding  the 
above  sources  of  non-differentiability  is  the  fact  that  many  of 
these  algorithms  use  upwind  algorithms,  which  have  their  own 
set  of  non-differentiabilities,  as  building  blocks.  In  particular, 
the  numerical  fluxes  associated  with  upwind  methods  are  non- 
differentiable  functions  of  the  hydrodynamic  wave  speeds,  at 
all  points  at  which  such  wave  speeds  vanish  in  the  frame  of 
reference  of  the  grid.  (These  are  the  places  in  the  flow  where 
u  —  c,  m,  and  u  +  c  vanish,  where  c  is  the  local  sound  speed 
and  u  is  the  local  fluid  speed.  We  shall  use  a  loose  defini¬ 
tion  here,  and  call  all  such  points  “sonic  points.”)  Thus  the 
use  of  even  simple  first  order  upwind  methods  like  Godunov’s 
method  [16]  would  seem  to  be  at  odds  with  our  desire  to  ac¬ 
curately  model  the  evolution  of  small  amplitude  perturbations, 
with  the  “high  resolution”  methods  being  riskier  still.  Let  us 
therefore  perform  some  numerical  experiments  to  see  if  we 
can  determine  just  how  risky  some  of  these  methods  may  be 
in  the  context  of  our  laser  pellet  implosion  problem. 

The  algorithms  we  choose  for  investigation,  in  order  of  in¬ 
creasing  differentiability,  are 

•  The  Piecewise  Linear  Method  (PLM)  of  Colella  and 
Glaz  [9]  (non-differentiable  whenever  the  slope  limiter 
is  active,  and  at  all  sonic  points) 

•  The  first-order  Godunov  method  [16]  (non- 
differentiable  at  all  sonic  points) 

•  A  first  order  donor  cell  method  plus  a  Richtmyer-Von 


Pellet  Material 

Pellet  Vapor 

Density  =  0.25  gm/cc 

Density  = 

0.0002  gm/cc 

Perturbed  Interface 


FIG.  4:  Geometry  of  laser-driven  pellet  test  problem. 

Neumann  artificial  viscosity  (non-differentiable  only  at 
entropy  wave  ( u  =  0)  sonic  points) 

•  The  Richtmyer  two-step  Lax-Wendroff  algorithm  ([17], 
pp.  302-303)(differentiable  everywhere) 

Note  that  only  the  Lax-Wendroff  method  is  truly  differen¬ 
tiable. 


V.  A  NUMERICAL  TEST  PROBLEM 

We  choose  as  our  test  case  a  simplified  version  of  the  early- 
time  pellet  implosion  problem.  As  depicted  in  Fig.  4,  a  semi¬ 
infinite  slab  of  a  monatomic  ideal  gas  at  one  atmosphere  pres¬ 
sure  and  a  density  of  0.25  g/cm3  is  irradiated  uniformly  from 
the  right  with  0.35  micron  laser  radiation,  with  an  intensity  of 
3  x  1013  W/cm2,  starting  at  t  =  0  and  lasting  indefinitely.  A 
simple  2D  code  that  we  dubbed  “Sandbox”  was  written  which 
incorporated  the  three  most  critical  physics  modules,  which 
were  applied  using  operator  splitting: 

•  Laser  absorption  via  inverse  Bremsstrahlung,  using  the 
prescription  of  Johnson  and  Dawson  [18],  modified  to 
ensure  differentiability. 

•  Spitzer-Harm  thermal  conductivity  [19]  without  flux 
limiting  (again  to  ensure  differentiability) 

•  one-dimensional  hydrodynamics  algorithms  applied  to 
two  dimensions  using  operator  splitting 

The  two-dimensional  calculation  was  then  performed  on  a 
800x16  grid  using  grid  spacings  of  0.1  microns  and  1.5  mi¬ 
crons  in  the  x  and  y  directions  respectively. 

The  interface  between  the  pellet  material  and  its  vapor  is 
placed  at  the  midpoint  in  a  cell,  ensuring  that  for  sufficiently 
small  interface  perturbations,  the  numerical  perturbation  will 
take  the  form  of  additive  density  perturbations  to  the  numeri¬ 
cal  grid  quantities.  (Recall  that  our  definition  of  differentiabil¬ 
ity  is  with  reference  to  additive  perturbations.)  A  sinusoidal 
perturbation  of  wavelength  24  microns  is  applied  to  the  po¬ 
sition  of  the  interface,  with  an  amplitude  of  0.0005  microns. 
This  amplitude  is  sufficiently  small  that  physical  nonlinear  ef¬ 
fects  should  be  extremely  small,  and  the  response  of  the  sys¬ 
tem  to  the  perturbation  should  be  linear  to  many  orders  of 
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magnitude.  The  physics  is  such  that,  after  some  initial  tran¬ 
sients,  a  quasi-steady  ablation  region  approximately  1.7  mi¬ 
crons  wide  establishes  itself,  propagating  into  the  pellet  ma¬ 
terial  as  ablated  low  density  material  is  blown  off  toward  the 
incoming  laser.  A  strong  shock  propagates  away  from  the  ab¬ 
lation  region  into  the  pellet.  Both  the  ablation  region  and,  to 
a  lesser  extent,  the  shock  are  regions  of  strong  perturbation 
evolution. 

Analysis  of  this  perturbation  problem  by  Velikovich  et  al. 
[20]  and  by  Goncharov  et  al.  [21]  reveals  that  the  evolu¬ 
tion  can  be  described  in  terms  of  a  quantity  that  can  be  mea¬ 
sured  experimentally,  the  integrated  mass  density  M(y )  along 
the  laser  path.  M(y )  starts  and  remains  a  constant  plus  a 
small  sinusoid.  Plotting  the  sinusoid’s  amplitude  as  a  func¬ 
tion  of  time  reveals  semi-periodic  nulls  at  which  phase  rever¬ 
sal  takes  place.  At  any  given  moment  in  time  during  a  calcu¬ 
lation,  we  should  be  able  to  compute  M(y),  and  compute  its 
Fourier  transform.  If  our  calculation  is  performing  properly, 
we  should  see  a  DC  component,  and  one  other  finite  amplitude 
component  at  the  seed  mode  wavelength.  The  amplitude  of 
all  other  Fourier  modes  should  be  due  solely  to  nonlinear  cou¬ 
pling  terms,  second  and  higher  order  in  the  perturbation  am¬ 
plitude,  and  hence  extremely  small  given  the  extremely  small 
perturbation  amplitudes  we  are  using  here.  Thus  in  the  fol¬ 
lowing  plots  we  use  the  following  definitions: 

•  “Seed  Amplitude”  is  simply  the  amplitude  of  the  seed 
mode  as  given  by  the  numerical  Fourier  transform  of 
M(y). 

•  “Noise”  is  the  result  of  numerically  removing  from 
M(y )  both  the  DC  component  and  the  seed  mode,  and 
taking  the  maximum  absolute  value  of  the  remaining 
values.  Note  that  this  is  an  extremely  sensitive  mea¬ 
sure  of  “noise.”  In  fact,  sufficiently  small  levels  of 
“noise”  may  in  fact  represent  physically  correct  second 
and  higher  order  nonlinear  coupling  terms,  as  we  men¬ 
tion  above. 


VI.  NUMERICAL  TESTS  OF  HYDRODYNAMICS 
ALGORITHMS  IN  2D  PLANAR  GEOMETRY 

The  Sandbox  code  has  been  used  to  run  many  tests  of  nu¬ 
merical  hydrodynamics  algorithms.  Here  we  have  chosen  just 
five  of  those  calculations  which  illustrate  the  potential  benefits 
of  maximizing  differentiability  in  such  algorithms. 

•  Test  1 :  Godunov  hydrodynamics  used  for  both  the  and 
y  directions. 

•  Test  2:  Godunov  hydrodynamics  used  for  both  the  x  and 
y  directions,  with  multidimensional  Colella-Woodward- 
Lapidus  artificial  dissipation  used  in  both  directions 
(see  discussion  below) 

•  Test  3:  Godunov  hydrodynamics  used  for  the  x  direc¬ 
tion  and  Lax-Wendroff  used  for  the  y  direction,  with  no 
additional  artificial  dissipation. 


FIG.  5:  Seed  mode  amplitude  (line)  and  noise  amplitude  (“+”  sym¬ 
bols)  versus  time  for  a  simulation  using  the  Godunov  algorithm  in 
both  the  x  and  y  directions. 

•  Test  4:  PLM  used  for  both  the  x  and  y  directions,  with 
no  additional  artificial  dissipation.. 

•  Test  5:  PLM  used  for  the  x  direction  and  Lax-Wendroff 
used  for  the  y  direction,  with  no  additional  artificial  dis¬ 
sipation. 

Test  1  In  Fig.  (5  we  show  the  seed  mode  amplitude  and 
noise  amplitude  versus  time  for  Test  1.  Although  the  seed 
mode  behaves  correctly  (see  below)  for  the  first  half  of  the 
simulation,  its  amplitude  is  thereafter  exceeded  by  that  of 
the  noise,  and  it  exhibits  incorrect  behavior  after  that  point. 
Our  diagnosis  is  that  the  Godunov’s  non-differentiability  at  all 
sonic  points  gives  rise  to  this  false  nonlinear  behavior,  making 
it  an  inappropriate  choice  for  the  accurate  modeling  of  small 
amplitude  perturbations  for  this  problem. 

Test  2  If  our  diagnosis  of  the  problems  seen  in  Test  1 
is  correct,  we  need  to  eliminate  some  of  regions  of  non¬ 
differentiability  in  our  algorithms  if  we  are  to  have  any  hope 
of  solving  our  problem  accurately.  We  will  do  that  below, 
but  first  let  us  look  at  another  possible  diagnosis  of  the  above 
problem:  that  we  are  witnessing  the  symptoms  of  the  “car¬ 
buncle  phenomenon”  [10,  23]  This  term  refers  to  the  diffi¬ 
culty  that  upwind  methods  like  the  Godunov  method  have 
when  dealing  with  perturbed  shock  waves  traveling  parallel 
or  nearly  parallel  to  one  of  the  mesh  lines  of  a  grid,  precisely 
the  situation  we  have  here.  While  the  true  causes  of  the  phe¬ 
nomenon!  are  still  a  subject  of  debate,  the  “cure”  is  nearly  al¬ 
ways  to  add  an  extra  amount  of  artificial  dissipation  transverse 
to  the  shock  propagation  direction.  Indeed,  Colella  and  Wood¬ 
ward  [10]  suggest  a  multidimensional  variant  of  Lapidus  arti¬ 
ficial  dissipation  which  closely  resembles  most  of  the  reme¬ 
dies  suggested  by  more  modern  papers.  Thus,  in  Fig.  6  we 
show  the  results  of  simply  adding  that  dissipation  to  the  algo¬ 
rithm  used  in  Test  1 .  We  see  that  the  “carbuncle”  hypothesis 
seems  to  be  plausible:  there  are  two  clean  phase  reversals,  and 
the  seed  mode  amplitude  is  decreasing  in  time,  in  accordance 
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FIG.  6:  Seed  mode  amplitude  (line)  versus  time  and  noise  ampli¬ 
tude  (“+”  symbols)  versus  time  for  a  simulation  using  the  Godunov 
algorithm  in  both  the  x  and  y  directions.  Extra  Colella-Woodward- 
Lapidus  artificial  dissipation  was  added  in  both  directions.  The  line 
labeled  “Peak  3  Amp”  indicates  the  expected  value  of  the  third  tem¬ 
poral  maximum  in  the  seed  mode  amplitude,  near  3  nsec,  using  the 
analytic  techniques  of  [22].  Note  that  although  our  results  show  good 
clean  linear  evolution,  with  the  expected  two  phase  reversals  in  the 
time  interval,  we  are  nonetheless  also  showing  considerable  false  nu¬ 
merical  damping  of  the  seed  mode  in  time. 

the  the  semi-analytic  results  of  [22],  which  predicts  a  slow  de¬ 
cay  of  the  seed  mode  in  time,  with  two  phase  reversals  in  the 
interval  between  0  and  3.0  nsec.  More  importantly,  except  in 
the  vicinity  of  the  phase  reversal  points,  our  signal-to-noise  ra¬ 
tio  is  roughly  105  to  1,  indicating  that  our  pure  Fourier  mode  is 
remaining  pure,  i.e.,  that  the  code  is  responding  in  the  differ¬ 
entiable/linear  manner  that  it  should.  However,  we  also  note 
in  the  figure  that  the  simulation  is  predicting  a  damping  of  the 
seed  mode  amplitude  in  time  that  is  considerably  larger  than 
we  expect  from  the  semi-analytic  results.  We  return  to  this 
point  later. 

Test  3  Instead  of  the  “carbuncle”  diagnosis  of  the  problems 
encountered  in  Test  1,  which  dictates  a  “cure”  of  increased 
numerical  dissipation,  let  us  now  pursue  a  solution  based  on 
our  hypothesis  of  non-differentiability  as  the  cause  of  Test  l’s 
difficulties:  We  know  that  the  only  non-differentiable  points 
for  the  Godunov  method  are  the  sonic  points,  which  exist  in 
both  the  x  and  y  directions.  We  also  know  that  our  perturba¬ 
tion  is  extremely  small,  meaning  that  the  “envelope”  in  so¬ 
lution  space  in  which  our  numerical  solution  will  encounter 
problems  will  be  correspondingly  small,  unless  of  course  the 
unperturbed  solution  itself  rests  directly  on  one  of  the  sonic 
points.  A  little  reflection  by  the  reader  will  hopefully  convince 
him  or  her  that  in  the  x  direction  we  would  have  to  be  quite 
unlucky  for  one  of  our  numerical  flux  evaluation  points  to  fall 
inside  the  envelope.  But  in  the  y  direction,  the  unperturbed 
solution  sits  directly  on  the  entropy  wave  sonic  point,  since 
for  the  unperturbed  solution  uy  =  0  everywhere.  Thus  if  our 
perturbed  solution  is  oscillatory  in  space,  which  it  is,  we  can¬ 
not  escape  being  caught  in  the  envelope  every  single  timestep. 


FIG.  7:  Seed  mode  amplitude  (line)  versus  time  and  noise  ampli¬ 
tude  (“+”  symbols)  versus  time  for  a  simulation  using  the  Godunov 
algorithm  in  the  x  direction  and  the  Lax-Wendroff  method  in  the  y 
direction.  No  extra  artificial  viscosity  was  added  in  either  direction. 
The  line  labeled  “Peak  3  Amp”  indicates  the  expected  value  of  the 
third  temporal  maximum  in  the  seed  mode  amplitude,  near  3  nsec, 
using  the  analytic  techniques  of  [22],  Note  that  we  have  achieved  the 
same  clean  results  as  in  Test  2,  and  at  the  same  time  have  been  able 
to  maintain  the  amplitude  of  the  seed  mode  near  the  value  given  by 
theory. 


This  reasoning  would  suggest  that  we  may  be  able  to  address 
Test  l’s  problems  by  using  a  differentiable  algorithm  in  the  y 
direction.  Since  the  solution  is,  for  these  early  times,  smooth 
in  that  direction,  we  can  and  do  choose  the  fully  differentiable 
Lax-Wendroff  method.  Note  that  we  are  suggesting  a  solu¬ 
tion  which  is  in  exactly  the  opposite  direction  from  that  given 
by  the  usual  “carbuncle”  remedy:  By  changing  from  the  Go¬ 
dunov  method  to  the  Lax-Wendroff  method  we  are  removing 
large  amounts  artificial  dissipation  rather  than  adding  it.  In 
Fig.  (7)  we  show  the  results  of  that  experiment.  Note  that  the 
results  are  even  better  than  those  of  Test  2.  Not  only  are  our 
signal-to-noise  ratios  roughly  10s  to  1,  but  the  seed  mode  am¬ 
plitude  near  3  nsec  is  much  higher,  in  very  close  agreement  to 
that  given  by  theory. 

We  again  point  out  that  we  have  arrived  at  two  totally  dis¬ 
parate  approaches  to  remedying  the  problems  of  the  Godunov 
method  in  Test  1:1)  adding  dissipation  in  the  y  direction;  and 
2)  removing  non-differentiability  in  the  y  direction,  and  in  the 
process  removing  dissipation  in  the  y  direction.  We  hypothe¬ 
size  here  that  the  first  remedy  is  nothing  more  than  a  cosmetic 
fix,  hiding  the  true  cause  of  the  “carbuncle”  phenomenon  un¬ 
der  a  mountain  of  numerical  dissipation.  Addressing  the  prob¬ 
lem  in  terms  of  its  true  cause,  i.e.,  failure  to  satisfy  the  differ¬ 
entiability  condition,  is  by  far  the  better  approach,  yielding  a 
clean  linear  response  without  the  side  effects  of  large  amounts 
of  numerical  dissipation. 

Tests  4  and  5  Although  we  do  not  include  the  figures  here, 
we  will  describe  the  results.  The  reader  may  see  that  Tests 
4  and  5  are  merely  the  higher  order  variants  of  Tests  1  and  3 
respectively.  They  yield  similar  results,  with  Test  4  showing 


unacceptable  growth  of  noise  in  the  time  interval  0  to  3  nsec, 
and  Test  5,  despite  the  added  but  improbably  crossed  non¬ 
differentiabilities  added  by  the  use  of  PLM  in  the  x-direction, 
yielding  good  clean  results,  albeit  with  signal-to-noise  ratios 
in  the  101 2 3 4  to  1  range,  as  opposed  to  the  105 6 7 8 9 10 11  to  1  range  of  Test 
3.  Thus  it  would  appear  that  even  the  more  modern  “high  res¬ 
olution”  methods  can  sometimes  be  used,  at  least  in  circum¬ 
stances  where  the  risk  incurred  by  their  failure  to  satisfy  the 
differentiability  condition  is  small.  The  user  of  such  methods 
must  still  be  wary  however,  for  they  can  and  will  occasionally 
trip  over  their  own  points  of  non-differentiability.  We  look 
forward  to  the  future  development  of  “high  resolution”  meth¬ 
ods  that  are  free  of  such  points.  Note  that  the  development 
of  “high  resolution”  methods  satisfying  the  differentiability 
condition  is  not  precluded  by  Godunov’s  Theorem,  since  the 
theorem  merely  demands  that  such  algorithms  be  nonlinear, 
not  that  they  be  non-differentiable. 

We  have  also  run  similar  tests  in  radial  geometry,  not  shown 
here,  wherein  shocks  can  focus  and  reflect  from  the  origin,  and 
we  find  there  an  increase  in  sensitivity  to  non-differentiability. 
We  have  found  at  least  one  case  where  our  perturbed  solution 
crosses,  albeit  “unluckily,”  a  u  +  c  sonic  point  in  the  radial 
direction,  causing  a  catastrophic  disruption  of  linearity  at  all 
subsequent  times.  However,  if  we  simply  reverse  the  sign  of 
the  1  =  2  Legendre  mode  perturbation,  we  get  “lucky”  again. 
We  alerted  the  reader  of  this  possibility  in  an  earlier  section. 
In  that  case  we  were  able  to  restore  proper  behavior  by  replac¬ 
ing  the  Godunov  algorithm  in  the  radial  direction,  which  has 
non-differentiabilities  at  all  sonic  points,  with  the  previously- 
described  donor  cell  algorithm,  which  is  non-differentiable 
only  at  u  =  0  sonic  points. 

VII.  CONCLUSIONS 

We  have  put  forth  the  hypothesis  that  a  highly  desirable 
condition  for  the  accurate  numerical  modeling  of  small  per¬ 


turbations  imposed  on  a  system  of  conservation  laws  is  that 
the  numerical  time  evolution  operator  satisfy  a  differentiabil¬ 
ity  condition.  The  numerical  experiments  that  we  present  here 
strongly  support  that  hypothesis.  Indeed  much  of  the  effort 
that  we  have  exerted  over  the  past  few  year  in  trying  to  ensure 
that  NRL’s  FAST  code  [3]  properly  models  small-amplitude 
perturbations  has  been  driven  by  this  design  principle.  We 
will  be  reporting  on  the  results  of  that  effort  elsewhere. 


However,  it  must  be  pointed  out  here  that  even  if  our  hy¬ 
pothesized  need  to  satisfy  a  differentiability  condition  is  cor¬ 
rect,  the  algorithms  we  have  used  here  to  meet  that  condi¬ 
tion,  i.e.,  using  a  different  hydrodynamics  algorithm  along  the 
shock  propagation  direction  than  transverse  to  it,  only  address 
the  needs  of  the  ICF  pellet  implosion  problem  in  the  early 
stages  of  its  evolution,  when  the  fronts  are  nearly  aligned  with 
the  transverse  direction.  As  the  perturbations  enter  the  fully 
nonlinear  regime,  the  need  for  a  front-capturing  method  in  the 
transverse  direction  will  probably  arise,  pushing  us  again  back 
in  the  direction  of  using  non-differentiable  algorithms.  Our 
empirical  experience  in  this  regime  suggests  that  the  situation 
is  not  as  bleak  as  it  sounds,  for  we  have  thus  far  found  differ¬ 
entiability  to  be  far  less  critical  an  issue  when  endeavoring  to 
model  the  evolution  of  large  amplitude  perturbations,  i.e.,  per¬ 
turbations  in  the  nonlinear  regime,  than  small  ones.  Nonethe¬ 
less  we  would  prefer  not  to  rely  on  this  experience.  Our  hope 
is  that  the  future  will  bring  forth  the  development  of  “high  res¬ 
olution”  methods  satisfying  the  differentiability  condition,  as 
we  discussed  in  the  last  section,  giving  us  the  robustness  of 
the  present  methods  in  the  presence  of  fronts,  combined  with 
the  differentiability  we  need  to  properly  model  the  evolution 
of  small  perturbations. 
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